%% CREATÉ TRIGGER TIME SERIES BY THRESHOLDING A TIME SERIES
% BUTTON PRESSES ARE TYPICAL APPLICATIONS BUT SOUNCARD OUTPUTS MIGHT BE USED TOO

%% CHECK WHICH ANALOGUE CHANNELS MIGHT BE INFORMATIVE
% Also get the time interval of interest
% ***User parameters
analogueChannels2Check={'ANIN1','ANIN2','ANIN3','ANIN4'};

%Visualize the analogue channels
ecogAna.data=[]
for k=1:length(analogueChannels2Check)
    load(analogueChannels2Check{k});
    ecogAna.data=[ecogAna.data; data];
end
baselineDurMs=0;
sampDur=1000/params.sampFreq;
ecogAna=ecogRaw2Ecog(ecogAna.data,baselineDurMs,sampDur,[]);
ecogTSGUI(ecogAna);

%% CREATE TRIGGER TIME SERRIES FOR ROM A SET
% ***User parameters
% Variables for trigger detection. Choice of thresholds might require some
% try and error
intervalStartEndSeconds=[16 230]; % defines the interval to use (in samples)
analogueChannels2Use={'ANIN3','ANIN4'};
minThresh=[0.1 0.1];; %An event is triggered when this value is eceeded. One entry for each analogue channel
minDurOfStimulusInSeconds=[0.5 0.5]; % That is the duration of the event.This is a "dead time" after an event was detected. Other events happening during this interval are ignored.

%variables for result visualization
checkResultVisually=1;
pre=0.5; % interval start prior to the event
dur=1; % displpay interval duration

%RUN THIS TO MAKE THE TRIGGERS

%Load the analogue channel
for k=1:length(analogueChannels2Use)
    
    load(analogueChannels2Use{k});
    intervalStartEndSample=nearly(intervalStartEndSeconds,(1:length(data))/params.sampFreq);
    data(1:intervalStartEndSample(1))=0;
    data(intervalStartEndSample(2):end)=0;
    %min trigger duration to samples
    minDurOfStimulusInSamples=round(minDurOfStimulusInSeconds(k)*params.sampFreq);
    threshData=double(data>minThresh(k));
    trigger=zeros(size(threshData));
    count=1;
    while count<intervalStartEndSample(2)
        if threshData(count)==1
            trigger(count)=1;
            count=count+minDurOfStimulusInSamples;
        else
            count=count+1;
        end
    end
    %Check the triggers
    %Actually we do this below after downsampling the trigger time series to
    %the sampling rate of the neuronal data
    if 0
        if checkResultVisually
            %check if all onsets are correct
            idx=find(trigger);
            figure;
            r=ceil(-pre*params.sampFreq:(-pre+dur)*params.sampFreq);
            timebase=r*1000/params.sampFreq;
            for m=1:length(idx)
                plot(timebase,[threshData(r+idx(m));data(r+idx(m));trigger(r+idx(m))*2]');
                input(num2str(m));
            end
        end
    end
    %Save what we have
    save(['trigger' analogueChannels2Use{k}], 'trigger', 'params')
    
    clear trigger params
    %Now we have to downsample to the same sample frequency as the brain data
    load Wav11
    wParams=params;
    load(['trigger'  analogueChannels2Use{k}])
    tParams=params;
    decFac=tParams.sampFreq/wParams.sampFreq;
    oldIdx=find(trigger>0);
    newIdx=round(oldIdx/decFac);
    data(:)=0;
    for m=1:length(newIdx)
        data(newIdx(m))=trigger(oldIdx(m));
    end
    
    params.sampFreq=wParams.sampFreq;
    save(['trigger'  analogueChannels2Use{k}],'trigger','data','params', 'tParams')
    
    
    %are the triggers good?
    load(analogueChannels2Use{k})
    [B,A] = butter(3,1/(2*decFac),'low'); %factor to nyquits of downsampled sequence
    data=filtfilt(B,A,double(data));
    aninDown=data(1:decFac:end);
    load(['trigger'  analogueChannels2Use{k}])
    save([ analogueChannels2Use{k} 'DownSamp'],'aninDown','params')
    
    load(['trigger' analogueChannels2Use{k}])
    % load the downssampled trigger and compared it to the downsampled analog channel
    idx=find(data);
    figure;
    %pre=0;dur=0.02;
    r=ceil(-pre*params.sampFreq:(-pre+dur)*params.sampFreq);
    timebase=r*1000/params.sampFreq;
    for m=1:length(idx)
        plot(timebase,[aninDown(r+idx(m));data(r+idx(m))*0.1]');
        title(['Channel:' analogueChannels2Use{k} ' trial # : ' num2str(m)])
        input(num2str(m));
    end
    if length(analogueChannels2Use)>k
        flag = 0;
        while flag==0
            r=input(['Continue processing next set? [y]/n '],'s');
            if strcmpi('y',r)
                flag = 1;
                % Next channel
            elseif strcmpi('n',r)
                error(['User aborted. Current channel: ' analogueChannels2Use{k}]);
                flag=0;
            else
                disp('Inavlid answer. Try again.')
                flag=0;
            end
        end
    end
    
end








